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We study the real-time dynamics of a two-dimensional Anderson-Hubbard model using nonequi¬ 
librium self-consistent perturbation theory within the second-Born approximation. When compared 
with exact diagonalization performed on small clusters, we demonstrate that for strong disorder this 
technique approaches the exact result on all available timescales, while for intermediate disorder, 
in the vicinity of the many-body localization transition, it produces quantitatively accurate results 
up to nontrivial times. Our method allows for the treatment of system sizes inaccessible by any 
numerically exact method and for the complete elimination of finite size effects for the times con¬ 
sidered. We show that for a sufficiently strong disorder the system becomes nonergodic, while for 
intermediate disorder strengths and for all accessible time scales transport in the system is strictly 
subdiffusive. We argue that these results are incompatible with a simple percolation picture, but are 
consistent with the heuristic random resistor network model where subdiffusion may be observed 
for long times until a crossover to diffusion occurs. The prediction of slow finite-time dynamics in 
a two-dimensional interacting and disordered system can be directly verified in future cold atoms 
experiments. 


Ergodicity plays a central role in the statistical me¬ 
chanics of closed systems . While it is usually difficult 
to rigorously prove that a given system is ergodic, our 
everyday experience strongly suggests that generic inter¬ 
acting systems with many-degrees of freedom are ergodic. 
The underlying assumption is that inelastic collisions be¬ 
tween particles allow for redistribution of energy and re¬ 
laxation to equilibrium. This is precisely the reason why 
it was commonly believed that Anderson localization [T], 
a truly nonergodic phenomenon, will be destroyed by the 
addition of a the smallest local inelastic interactions [5]. 
It therefore came as a surprise when two groups argued 
that a truly nonergodic phase, later dubbed the many- 
body localized phase, generically exists at a finite energy 
density and finite interaction strengths [111]. By tuning 
the parameters of the system a dynamical many-body 
localization transition between ergodic and nonergodic 
phases occurs. This transition was first observed numer¬ 
ically EllS], and more recently in cold-atoms experiment 

El- 

While the existence of a stable nonergodic phase at 
finite energy density is the hallmark of many-body local¬ 
ization, transport, or the absence of thereof, in this phase 
is not very interesting. On the other hand, it was recently 
demonstrated that transport within the ergodic phase, 
previously believed to be diffusive El SI, exhibits subd- 
iffusion for nontrivial times in one-dimensional systems 
EHin]- A heuristic theoretical explanation of the subd- 
iffusive behavior suggests that transport in the system 
is similar to the random resistor problem [111 I12j , where 
rare and large insulating regions have a dominant effect in 
one-dimension m- Subdiffusion also naturally arises in 
the phenomenological renormalization group (RG) treat¬ 
ments of the same problem, where close to the MBL 
transition blocking (or thermal) regions appear on every 
length scale (RG step) [131114] . The fractal structure of 
the blocking (or thermal) regions has led to the proposal 


that the many-body localization transition is a sort of a 
percolation transition |14til6| . Similar suggestions were 
previously made for the metal-insulator transition in two 
dimensional systems at zero temperature ini- 

One of the predictions within the percolation picture 
is that the observed subdiffusive behavior in the ergodic 
phase in one-dimensional systems will be absent in higher 
dimensions, where blocking regions can be avoided |14L 
in EH]. This prediction is asymptotic in its nature, and 
any generalization to finite times will depend on the cho¬ 
sen model. One of the main problems in studying subd¬ 
iffusion is that most numerically exact methods struggle 
within the ergodic phase. This is in contrast to the case 
of the nonergodic MBL phase, where many numerically 
exact methods are available and efficient (see |191 ED] and 
references within), . One exception is the study of dc con¬ 
ductivity at low temperatures using determinantal quan¬ 
tum Monte-Garlo (DQMG), where it was possible to cir¬ 
cumvent the need for analytic continuation and compute 
dc conductivity directly in real time. This allowed for 
the study of the metal-insulator transition as a function 
of the interaction strength for a two-dimensional model 
|21H23| . It is not clear, however, if this approach carries 
over to high temperatures, of interest to MBL. Moreover, 
it provides only asymptotic (dc) information. 

In this work, we examine the dynamics of a two- 
dimensional Anderson-Hubbard model using nonequilib¬ 
rium perturbation theory within the second-Born ap¬ 
proximation. While our approach is approximate, we 
show that it quantitatively reproduces exact results in 
small systems for a wide range of parameters and ap¬ 
pears to become exact for the timescales of investigation 
in the limit of strong disorder. This approach is similar to 
the diagrammatic approach of Ref. |3], which first estab¬ 
lished the many-body localization transition. We relax 
several of the approximations of Ref. [^ and compute in 
detail the dynamics of the system from an appropriately 


2 


chosen initial condition. 

We investigate a two-dimensional Anderson-Hnbbard 
model, 


H — t ^ ^ “1" U ^ ^ d" ^ ^ ^rcrhrCT^l) 

(rr'),<T r r.cr 

where creates a fermion of spin a = ±1/2 at site r, 
hrcr is the density operator, t is the hopping matrix el¬ 
ement (we set t = 1 thronghout), U is the interaction 
strength and hie, are random fields independently dis¬ 
tributed on the interval hia G [—W, W]. To eliminate all 
symmetries of this Hamiltonian we subject the fermionic 
species to different disorder fields. Since the model has 
a finite number of states per site the energy density is 
bounded, which allows for meaningful consideration of 
the infinite temperature limit. Note that while for sys¬ 
tems with unbounded energy density the infinite temper¬ 
ature limit normally coincides with the classical limit, 
this is not necessarily the case for systems with bounded 
energy density. For example in the noninteracting Ander¬ 
son model in one and two dimensions, all the eigenstates 
are localized and therefore localization persists at any 
temperature. In this limit the many-body localization 
transition occurs as a function of the parameters of the 
system, and not the temperature [S]. In what follows we 
will consider only the infinite temperature limit, and use 
the strength of the disorder as a control parameter of the 
many-body localization transition (we fix the interaction 
to be f7 = 0.5). The theoretical boundary of the MBL 
transition at half-filling and infinite temperature may be 
estimated as (Eq. (93) of Ref. |24|'). 

= K = 4A.C^ (2) 

where v, A and ^ are the one-particle density of 
states, bandwidth and localization length, correspond¬ 
ingly. Taking the appropriate values for the one-particle 
problem, the boundary line of the MBL transition can be 
readily calculated (see Fig. [^. Since we are interested in 
studying the dynamics across the transition, we vary the 
disorder strength in the interval 5 < kF < 20 (yielding a 
non-interacting localization length ^ < 5.5). 

To study the dynamics we use the equations of motion 
for the one-particle nonequilibrium correlators, 

Gtj (^; t') = -*Tr |/5oc* (t) c\ (t')} , (3) 

Gfj [t] t') = iTr jpoc] {t') Ci (t)| , 

where po is the initial density matrix. For an uncorre¬ 
lated initial density matrix, the Green’s functions obey 



Figure 1: (color online) Dynamical phase boundaries of the 
two-dimensional disordered Hubbard model at infinite tem¬ 
perature as a function of W and U . The phase boundary was 
obtained using Eq. (93) of Ref. [24] and demarcates ergodic 
and nonergodic regions. The white circles correspond to the 
parameters used in this work IT = 5 — 20 and U = 0.5. 

the Kadanoff-Baym equations of motion [25], 

idtG^itff') = + {t))G^{tff') 

± f i:^{t,t2)G^ihff')dt2 

Jo 

nt' 

+ / T.^{t,t2)G^{t2,t')dt2, (4) 

Jo 

where spatial indices and summations are suppressed 
for clarity; /ig is the one particle Hamiltonian; (t), 
(t) are the Hartree-Fock greater and lesser self¬ 
energies of the problem respectively; and the superscripts 
’R’ and ’A’ represent retarded and advanced Green’s 
functions and self-energies, which are defined as 

E«(t,t2) = 0(<-t2)(S>(t,<2)-S<(t,t2)) (5) 

G^(f2,t') = -e{t'-t2){G>it2,t')-G<{t2,t')). 

We note in passing that this equation is exact if the exact 
self energies are employed. Since normally this is not 
possible, approximate forms for the self energies are used. 
Here we use the second-Born approximation for the self¬ 
energy, 

= -iU5,,G<{t-t) 

E>(t,t') = U^G<{t\t)G>{t,t')G>{t,t') (6) 

which is a self-consistent conserving approximation, 
specifically it conserves the total energy and number of 
particles and amounts to a resummation of an infinite 
class of diagrams. This approximation was originally 
considered in the pioneering work of Ref. jl] , but due to 
the complexity of it was reduced to a corresponding 
quantum Boltzmann equation. It is numerically feasible 





3 


to solve Q, within full second-Born approximation, how¬ 
ever this is known to generate spurious relaxation |26| . 
This problem can be avoided via the the introduction 
of the generalized Kadanoff-Baym anzatz for the greater 
(lesser) correlator, 




(t, t') {f, t') - G^ (t, t) G^ (t, t') , 


(7) 

and by approximating the retarded and advanced Green’s 
functions with their Hartree-Fock (HF) values |27fl29| . 
Substituting Eqs. Q and ^ into Q one obtains the 
Quantum Master Equation (QME) for the one-particle 
density matrix.Unlike the Boltzmann equation used in 
Ref. |3] , this approach considers the full one-particle den¬ 
sity matrix and not just its diagonal values, and therefore 
includes additional quantum information. We refer the 
reader to Refs. p7l [28] for additional discussion. 

In the infinite temperature limit, the equilibrium den¬ 
sity matrix po iii is proportional to the unity matrix, 
and it is tempting to use this maximally mixed state as 
an initial condition. However as alluded to in Ref. U and 
we explained in our previous work jHj, this choice would 
lead to erroneous results, in particular in the limit of in¬ 
finite disorder. To alleviate this issue we use pure and 
uncorrelated initial conditions. Specifically, we fix the 
number of fermions for each species (we use half-filing 
throughout this work) and distribute all fermions ran¬ 
domly on a square lattice. For these initial conditions 
it is straightforward to show that our approach becomes 
exact in the zero hopping limit. Moreover, as we demon¬ 
strate below, even for finite hopping and large disorder 
our method quantitatively approaches the exact behavior 
exhibited by finite sized systems studied with exact di- 
agonalization. The infinite temperature limit is achieved 
by averaging the observables over random realizations of 
the initial conditions. To study the dynamics of the sys¬ 
tem we calculate the density-density correlation function 
at infinite temperature. 


C{r,r';t) 



( 8 ) 

where |a) states are the random initial states described 
above and Z is the total number of many-body states. 
Normally this quantity requires two-particle correlators, 
however our initial conditions have the useful property, 
(h'ro' (f) ^r'(7 ( 0 )) — (h-rcT (f)) ^r'tr (b)) where ( 0 ) — 
0,1 is the occupation number of the initial state and. 


{fira (t)) = -iG< (rt, rt). (9) 


To study the relaxation of the system we calculate the 
autocorrelation of the density and the mean-square dis¬ 
placement, 

r r 

(10) 



Figure 2: (color online) Comparison between exact diagonal- 
ization (black) and nonequilibrium self-consistent perturba¬ 
tion theory (red/gray) for several disorder strengths {W = 1,5 
and 7) performed on a small Hubbard cluster of dimensions 
3x3. The left panel shows the density-density correlation 
function as a function of time, and the right shows tlm aver¬ 
age mean-square displacement, computed using Eq. (lOl. An 
average over 1000 disorder realizations was performed and 
shaded areas designate uncertainty bounds. 


For a diffusive system we have, p {t) ^ and (t) ~ 

t, while for localized systems both p (t) and (t) are 
expected to saturate at infinite times to finite quantities. 
Similar quantities based on the one-time density [rir (t)] 
were studied for the clean systems out-of-equilibrium isoi- 
I32| . however, such quantities cannot be directly used in 
equilibrium where any one-time operator is independent 
of time. 

The quantum master equation (QME) is solved nu¬ 
merically for the one-particle correlators as a function of 
time within the self-consistent second-Born approxima¬ 
tion. For this purpose we utilize the parallel numerical 
scheme developed in Ref. |33]- After obtaining the one- 
particle correlators, we extract the density-density cor¬ 
relation function and average it over both the initial 
j conditions and disorder realizations. Depending on the 
strength of the disorder, we use between 35 realizations 
for the weakest disorder strengths and 400 realization for 
the strongest disorder strengths. The perturbation the¬ 
ory reproduces the exact (and trivial) result for zero hop¬ 
ping. To evaluate its performance for nonzero hopping 
we compare the perturbative calculation to the exact di- 
agonalization (ED) solution of a small Hubbard cluster 
of size 3x3. As seen from Fig. a remarkable corre¬ 
spondence between the QME and the exact solutions is 
seen for U = 0.5 and times t < AO. This correspondence 
does not exist at the HF level (not shown). While this 
comparison has little bearing on the dynamics of MBL 
in the thermodynamic limit we will assume that this cor¬ 
respondence does not become worse for larger systems 
sizes. 

The nature of the transport is assessed by examining 
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the finite time dynamical exponent, 


a{t) = 


dlnr^ (t) 
dint 


( 11 ) 


which has values a (t —>■ oo) = 2 for ballistic transport 
and a (t —>■ oo) = 1 for diffusive transport. This exponent 
is directly related to the dynamical exponent, z = 2/a, 
which is commonly used in the context of phase transi¬ 
tions and defines the scaling between time and space, 
t (X . For finite systems, asymptotic time dynam¬ 
ics will be dictated by finite size effects such as reflec¬ 
tions from the boundaries and the asymptotic value of 
a (t —> oo) will have little meaning. Nevertheless, due 
to the locality of the interaction it takes a finite time, 

, for a local perturbation to reach the boundary of the 
system [Mj. Therefore for times t finite-size ef¬ 

fects can be effectively eliminated. In this work we have 
chosen = 50, which is sufficient to remove finite size 
effects for the chosen parameters and sizes studied. Such 
calculations required about 300, 000 computer hours For 
the lowest disorder strength we investigated systems of 
size 24 X 24 were required while for the slowest dynamics 
(strongest disorder strength) systems of the size 16 x 16 
were sufficient. The reader is referred to our previous 
works for a more detailed explanation of this procedure 
(HIS]. We note in passing that the required system sizes 
lie far above the accessible sizes of any exact numeri¬ 
cal method. Our method is approximate, however Fig. 
suggests that for the studied disorder strengths and con¬ 
sidered times it is essentially quantitatively exact. While 
the chosen time scale of t* = 50 is definitely not asymp¬ 
totic, as we will show below it allows access to nontrivial 
dynamics across the MBL transition. Moreover this is 
the relevant time-scale for current experiments in cold- 
atoms, where finite size effects from the harmonic trap, 
particle loss and decoherence set-in on longer time scales 
[ 21135 ]. 

To evaluate the dynamical exponent as a function of 
time [see 01 , we divided the data into equal time in¬ 
tervals (on a logarithmic scale) and performed piecewise 
linear fitting. This is a well defined procedure, which 
produces a time-dependent dynamical exponent, and is 
different from numerical extrapolation performed using 
power-law fitting. Both procedures however have little 
bearing on asymptotic transport. The result is demon¬ 
strated for several disorder values in the left panel of 
Fig-i For localized systems the value of the dynami¬ 
cal exponent as a function of time, a (t), should decay 
to zero, while for diffusive systems it should asymptot¬ 
ically saturate to a (< —>■ oo) = 1. As clearly seen from 
Fig-i neither occurs for the studied system, instead the 
dynamical exponent saturates to a finite value on the 
studied timescale. We repeat this analysis for various 
disorder strengths and show that the dynamical expo¬ 
nent converges to a finite plateau value for all cases un¬ 
der investigation in the right panel of Fig. i This fig¬ 



Figure 3: (color online) Dynamical exponent across the many- 
body localization transition. On the left panel the dynamical 
exponent is determined as a function of time for different dis¬ 
order strengths. Darker shades designate stronger disorder. 
The right panel shows the dynamical exponent as a function of 
the disorder strength. Darker shades represent longer times, 
and the arrows illustrate the convergence of the dynamical 
exponent in time. 


ure summarizes the main result of our work: for non¬ 
trivial times the dynamical exponent is finite across the 
many-body localization transition and deep in the er- 
godic phase. This means that for experimentally rele¬ 
vant times and a broad range of parameters the system 
exhibits subdiffusive transport. While the transition to 
diffusion can and may occur at later times (unaccessible 
to us) , as happens in the case of classical supercooled 
liquids |3S], we do not observe this crossover in our study 
even for the lowest disorder strength [W = 5]. Studying 
even lower disorder strengths is problematic due to the 
exponential dependence of the non-interacting localiza¬ 
tion length on the disorder strength in two-dimensional 
systems, ^ oc exp [A/VF^] |3Zj. From our data it is diffi¬ 
cult to pinpoint the location of the MBL transition since 
a ft) is never exactly zero. Nevertheless for W = 15 — 20, 
the values of the dynamical exponent we obtain are zero 
within the error bounds, which is consistent with the 
theoretical estimate presented in Fig. We empha¬ 
size that our results do not contradict the expectation 
that the MBL dynamics within the ergodic phase of one¬ 
dimensional systems is qualitatively distinct from that 
in higher dimensions, however the duration of apparent 
subdiffusion we observe is surprising. 

Note that our results are rather different from dynam¬ 
ics close to the noninteracting Anderson transition |38j . 
In this case, close to the transition, the system exhibits 
apparent subdiffusion with a constant dynamical expo¬ 
nent a = ^jd, for lengths of the order of the correlation 
(localization) lengths and then crosses-over to either dif¬ 
fusion (a = 2) or localization (a = 0). This hints that 
the MBL transition is inherently different from the An¬ 
derson transition. It is also interesting to understand 
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how our results compare to the heuristic percolation pic¬ 
ture USHIS]- For classical percolation, the dynamics of 
the system at the percolation threshold is subdiffusive 
on all time scales with the universal dynamical exponent 
\a « 0.69, for d — 2], Above the percolation threshold 
a (t) monotonically increases until it reaches one (nor¬ 
mal diffusion), and below the percolation threshold a (f) 
monotonically decreases to zero |39]. Our results there¬ 
fore can only correspond to the nonergodic phase below 
the percolation threshold. This renders a literal mapping 
of the MBL transition to a classical percolation transi¬ 
tion as suggested in Refs. usmn inconsistent with the 
theory of Refs. 110121 ], which places the W = 5 point 
deep within the ergodic phase (see Fig. [^. On the other 
hand the random-resistor network picture cni can be 
consistent with our results. Within this picture for a 
sparse distribution of the resistors, an initial monotonic 
decrease of a (t) can display long plateaus and an even¬ 
tual crossover to diffusion |40]. While we have not ob¬ 
served this crossover even for the lowest studied disorder, 
we do observe long plateaus of a (t) which can signal a 
precursor of this behavior. 

In summary we have studied the nonequilibrium dy¬ 
namics of a two-dimensional Anderson-Hubbard model. 
For this purpose we have utilized self-consistent nonequi¬ 
librium perturbation theory. Our method allows for the 
study of large system sizes and the elimination of finite 
size effects up to nontrivial times. While our method 
is approximate, it appears to become exact in the limit 
of strong disorder and gives quantitatively reliable re¬ 
sults even for intermediate disorder strengths when com¬ 
pared to exact diagonalization studies on small systems. 
We present evidence consistent with the existence of the 
many-body localized phase at a two-dimensional system 
at infinite temperature, and put an upper bound on the 
location of the many-body localization for a particular in¬ 
teraction strength. Surprisingly, up to the studied times 
we find subdiffusive dynamics for all the studied param¬ 
eters. We compare our results to the predictions of the 
available heuristic theories and point out that while they 
can be explained within the random-resistor network pic¬ 
ture they are inconsistent with the literal adoption of 
classical percolation theory for the MBL transition. Our 
results predict that experiments performed on cold atoms 
in two-dimensional disordered systems will observe a fi¬ 
nite regime of subdiffusive transport for times shorter 
than the typical decoherence time in these experiments. 

We would like to thank G. Cohen for many enlighten¬ 
ing and helpful discussions. This work was supported by 
grant NSF-CHE-1644802. 
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